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ABSTRACT 

y  To  simulate  long-run  averages  of  time  integrals  of  a  recurrent  semi- 
Markov  process  efficiently,  convert  to  discrete-time  by  simulating  only  an 
imbedded  chain  and  computing  the  conditional  expectations  of  everything  else 
needed  given  the  sequence  of  states  visited.  This  reduces  asymptotic  variance 


and  eliminates  generating  holding-time  variates.  In  this  setting, 

77 1 

uniformizing  continuous-time  Markov  chains  is  not  worthwhile.  -We- generalize 


beyond  semi-Markov  processes  and  cut  ties  to  regenerative  simulation 


methodology. 
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SIGNIFICANCE  AND  EXPLANATION 


A  broad  class  of  stochastic  systems  which  are  studied  in  operations 
research  may  be  modelled  as  semi-Markov  processes.  Frequently,  one  is 
interested  in  obtaining  an  estimate,  via  simulation,  for  the  steady-state 
average  of  such  a  process.  In  this  paper,  we  offer  new  insights  on  an  easily 
implemented  procedure,  which  can  substantially  improve  the  accuracy  of  such  an 
estimate.  The  basic  idea  involves  passing  from  the  continuous-time  semi- 
Markov  process  to  an  appropriate  discrete-time  sequence,  by  conditioning  out 
holding  time  variables. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 
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DISCRETE-TIME  CONVERSION  FOR  SIMULATING 
SEMI-MARKOV  PROCESSES 

13  2  4 

Bennett  L.  Fox  '  and  Peter  W.  Glynn  ' 

1 .  Introduction 

Let  X  be  a  positive  recurrent  semi-Markov  processes  and  f  be  a 
real-valued  function.  We  want  to  estimate  the  time  average 

r  -  11m  ( 1/T) E  JT  f(x(t))dt  (1) 

0 

assuming  It  exists.  Section  2.1  starts  with  a  standard  regenerative 
simulation  approach  and  then  converts  to  discrete  time  by  conditioning  on  the 
sequence  Y  ■  (Y0,Yj  ,...)  of  states  visited  in  an  imbedded  Markov  chain.  This 
reduces  variance  and  also,  typically,  the  work  to  simulate.  Hordijk, 
Iglehart,  and  Schassberger  [6]  adopt  the  same  approach  for  the  special  case  of 
continuous-time  Markov  chains,  except  that  they  prove  variance  reduction  by 
explicit  calculation  without  mentioning  the  general  principle  that  computing 
conditional  expectations  reduces  variance.  Peter  Lewis  informed  us  that  he 
too  was  aware  that  this  conditional  Monte  Carlo  approach  would  streamline  the 
proofs  in  [6].  Fox  and  Glynn  [5]  obtain  an  analog  of  these  results  for 
certain  finite-horizon  semi-Markov  processes.  Unlformizatlon  pays  in  [5]  but 
not  here  as  shown  in  [6]  and,  with  less  effort,  in  section  3. 

Section  2.2  considers  more  general  processes  and  cuts  ties  to  the 
regenerative  approach. 
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We  reach  state  Yj  at  tine  Tj.  Let 


«.  ■  1,  -  T, 

k  k+1  k 


(10) 


P(ak  e  dt  |  Y)  -  F(Yk,Yk+1,dt) 


y(x,y)  -  J"  zF(x,y ,dz) . 
0 


(11) 

(12) 


As  jump  Nr  ends,  we  visit  state  i  for  the  n-th  time.  This  gives 


N  -1 
n 


?„-(*/»)  io  .-(vv.) 


(13) 


N  -1 
n 


E[Vn  |  Y]  -  (  1/n)  ^  fCYjJwCYj.Y^j)  , 


(14) 


so  we  can  compute  r  in  (9).  Computing  and  accessing  the  expected  transition 


times  and  the  transition  probabilities  are  similar  tasks.  Fox  and 

3 


Glynn  [5]  and  references  cited  there  discuss  the  latter.  When  p(y.,Y.^,) 


depends  only  on  Y^,  the  former  job  is  normally  easy. 


Since  X  is  a  semi-Markov  process  (by  assumption),  the  <»k' s  are 
conditionally  Independent  given  Y.  So  Y  regenerates  implies  X  regenerates. 
Assuming  certain  mild  moment  conditions  [Inequality  (24)  will  do]  and 
mimicking  standard  arguments,  we  get 

4(rn  -  r)  ->  (6/EtJ  N(0,1)  (15) 

/»(rn  -  r)  ->  (o/EtJ  N(0,1)  (16) 

where 


3 


>  V 


62  -  Var  07) 

o2  -  Var  (e[d^  I  y]}.  (18) 

Using  the  standard  variance-reducing  property  of  conditional 

expectation,  we  get  a  <  6.  Clearly,  Etj  -  ETj  ,  so  rn  has  smaller  variance 

A 

than  r  . 
n 


2.2  Generalization  and  formal  proof 

Now  we  neither  assume  that  Y  is  a  Markov  chain  nor  that  the  a  are 
conditionally  Independent  given  Y.  The  only  structural  assumptions  relating  Y 
and  X  are  (11)  and  (19).  With  these  provisos,  define  Yfe,  a^,  and  y(x,y)  as 
before.  Let  1  be  an  indicator,  TQ  -  0,  +  •,  and 

00 

X(t)  -  l  Y  l(T  <  t  <  T  ).  (19) 

k-0  * 

Thus ,  X  is  more  general  than  a  semi-Markov  process. 


The  obvious  estimator  is 
T 

R  -  (l/T  )  /  n  f(  X(s))ds  (20) 

n  n  0 

Typically,  the  expected  work  to  generate  Rr  or  rn  grows  at  rate  n  whether  n 
indexes  "transitions"  as  in  this  section  or  cycles  as  in  section  2.1.  So  we 
compare  efficiencies  of  various  estimators  with  respect  to  n. 


Roughly  speaking,  to  convert  to  discrete  time  we  compute  conditional 
expectations  given  Y.  The  alternative  estimator 
n-1  n-1 


■  l0 


) 


(21) 


is  certainly  plausible.  The  only  difference  between  Rr  defined  by  (20)  and  rQ 


defined  by  (7)  lies  in  what  n  Indexes;  likewise ,  for  Rq  and  rQ.  We  show  that 


-srfVfr. 


■< '  •*  *-  * 


Rn  beats  Rm  in  a  precise  sense,  provided  that  Rr  Is  no  harder  to  compute  than 


We  assume  that  X  satisfies 


(A)  T  /n  ->  B  >  0 
n 


(6  constant) 


(B)  there  are  finite  constants  rt  and  Oj  such  that 

(1/4)  /  n  f(x(s)-ri)ds  ->  a1  N(0, 1) 

0 

(C)  there  are  finite  constants  r^  and  o2  such  that 

r  n-1 

(l^n)  l  (f(Tk)-r2)p(^kiYk+1)  «>  o2  N(0,1) 

-1 

(D)  the  sequences  { n  £  f(Y.)o.:  n  >  l}  , 

k-0  *  * 

T 

{ T _/n:  n  >  l}  ,  and  {  n  1  [/  n  (  f(  X(s))  -r)  ds]2  :  n  >  l} 

0 

are  uniformly  integrable. 


Section  2.3  contains  examples. 


Not  surprisingly,  we  have 


Proposition.  Tj  “  r2  • 


i  T 

1  r  n. 


Proof.  From  (A)  and  (B),  we  get  ~J  nf(x(s))ds->  r  B.  Equivalently, 
t  n-1  0  1 

—  I  f(Y,  )a  ->  r,B.  From  (D),  we  also  get  weak  convergence  to  r  a  in  the 
k=0  k  x  i  l 

function  space  Lj  .  From  p.  306  of  Chung  [4],  we  therefore  get 
1  0-1 

-  I  E( f ( Y  ) a  I  Y)  ->  r  B; 
k-0  K  * 

*•••  n  °L  ->  ri6- 


mw- 


The  left  side  is  the  numerator  of  R  .  Likewise,  for  the  denominator  of  R  we 

n  n 


*2  u(vYk+i>  ->6- 

From  (C),  R^  *>  r2 .  Comparing  these  results  proves  that  rj  ■  r2  .  I 


We  are  now  ready  to  compare  R  and  R  .  From  (A)  and  (B) : 

n  n 


4  (Rn-ri)  «>  (Oj/B)  N(0, 1)  ; 


from  (A)  and  (C): 


4  (R  -r  )  ->  (o  /B)  N(0 , 1 )  , 


That  confidence  intervals  with  fixed  coverage  based  on  R  are 

n 

asymptotically  shorter  than  those  based  on  R^  now  follows  from  our  main 
2  2 

Theorem.  o2  <  Oj  . 

Proof .  By  Jensen's  inequality  for  conditional  expectations  (see  p.  302  of 
Chung  [4]), 


i  T 

r  nr 


"  'e(|  n(  f(  X(s)) -r)  ds)2  >  n  1E[(/  n(  f(  X(s)) -r)  ds  |  y)Z] 

0  0 

“  n_1  e[  ^  (ffYjJ-rJuCYj.Y^)]2. 

2 

By  (D)  the  extreme  left  side  converges  to  c  .  We  show  that  the  right  side 
2 

converges  to  ,  which  requires  a  uniform  integrability  argument. 


Av  "  l""1!  I  (f(Yj-r)w(Y4,Y^,)]2  >  k}. 


Since  A^  is  Y-measurable ,  apply  Jensen  s  inequality  again  to  get 

*UA  («(Y.)-r)*(Y..Twl))2}  <  E[(I.  (i/T”  f(x(s)-r)ds)f]. 

nk  j-0  J  J  nk  0 

Use  (C)  to  see  that  lim  sup  p(  A  0  and  then  theorem  4.5.3  of  Chung  [4]  to 

k+®  n 

see  that  the  term  in  brackets  goes  to  0  uniformly  in  n  as  k  +  ®.  This  proves 
uniform  integrability  of  the  term  in  braces.  II 


2.3  Examples 


1.  Let  X  be  a  semi-Markov  process  for  which  Y  is  regenerative:  there  exists  s 

such  that  P[Y  *=  s  infinitely  often]  *  1.  Set  T  ■  inf{  n  >1:  Y  =  s}  and 
n  n 

assume  that  YQ  ■  s  and 
T-l 

E(  I  (|f(Y)|  +  l)a)  <».  (24) 

k=0  K  K 

Then  (A)-(D)  hold;  for  uniform  integrability,  see  Chung  [3]. 

2.  Let  X  be  a  semi-Markov  process  for  which  Y  is  a  stationary  ♦  -mixing 
process.  If  Y  is  ♦-mixing,  then  l^(Yn)an}  is  ♦-mixing  with  mixing 
coefficients  doubled  since 

|P{f(Y0)a0  e  A,  f(Yn)<»n  e  »}  -p{f(Y0)a0  e  A} -p{  f(  yJ^  e  b}  | 

-  |E{p{f(Y0)aQ  e  A  |  Y}*p{f(Yn)an  c  B  |  y}  }  -  p{  f(  YQ  )°to  e  A}  •  p{  f(  yJ  e  b} 

-  -%(VY,)-'Mw1)' 

(where  gj(x,y)  -  p{  f  (  TCq  )  <»0  e  A  I  Y(  ■  x,  T(  -  y} 


gjCx.y)  -  P(ftvjan  e  B  |  -  x,  Yn+[  -  y) ) 


‘  2*nEt,(’'o-Yl)  '  2VP(f(Y«H  e  *!• 


’  **  *«  • 


PtlTf 


■  y~s\  ■' .» w 


MrtrCrCifiSWl 


the  inequality  following  from  (20.28)  of  Billingsley  [1].  Suppose  that 


and  E[(f(Y0)+l)a0]2  <«. 

Then,  (A)  -  (D)  hold  (see  theorem  20.1  of  Billingsley  [1]). 


3.  Unifonaization  loses 


For  the  special  case  of  continuous-time  Markov  chains  with  conservative 
generator  Q  and  jump  rates  having  least  upper  bound  X  <  •  ,  we  can  unlformize. 
This  corresponds  to  the  equal  holding-time  method  of  [6].  The  method  of 
section  2.1,  without  unlforml zation ,  corresponds  to  the  constant  holding-time 
method  of  [6].  With  uniformization ,  the  imbedded  discrete-time  Markov  chain  W 
has  null  jumps  from  a  state  to  itself.  Thus,  to  generate  any  fixed  number  of 
regeneration  cycles  with  the  equal  holding-time  method  requires  more  jumps, 
hence  more  work.  It  also  gives  more  variance,  as  shown  in  [6]. 

We  now  give  a  simpler  proof.  Delete  all  null  jumps  from  W  to  get  Y, 
giving  o(Y)  co(W)  where  o(H)  is  the  o-field  generated  by  H.  By  proposition 
G.l  in  [5]  for  example,  we  get  higher  variance  by  conditioning  on  the  latter. 
Thus 

Var  (Efl^  |  y]  )  <  Var  (  e[  1^  |  w] ) .  (25) 

Nevertheless,  what  is  the  best  way  to  uniformize?  To  get  a  legitimate 
representation,  we  must  choose  the  (Poisson)  clock  intensity  0  >  X  in  the 
uniformized  process.  Choosing  0  ■  X  stochastically  minimizes  the  minber  of 
jumps  in  W  to  simulate  a  fixed  number  of  regeneration  cycles,  hence 


stochastically  minimizes  work.  Hordijk,  Iglehart,  and  Schassberger  [6]  show 

2 

by  explicit  calculation  that  choosing  6  •  X  also  minimizes  o  in  (18).  An 
outline  of  an  easier  proof  follows.  Subscript  W  and  X  to  indicate  clock 
intensity.  Put  X^  and  X0  on  the  same  probability  space  using  three 
synchronized  streams  of  common  random  numbers.  Stream  1  generates  the  common 
non-null  jump  subsequences  of  W.  and  W.  .  Stream  2  generates  the  null  jumps  so 
that  each  null-jump  sequence  in  is  at  most  as  long  as  the  corresponding 
null-jump  sequence  in  W0 .  Stream  3  generates  clock  chimes  so  that  chime  j  of 
the  0-clock  sounds  before  or  with  chime  j  of  the  X-clock.  Appendix  B  of  [5], 
in  a  more  general  setting,  spells  out  the  details.  Thus 

o(wA)  c  o(w0)  (26) 

for  all  0  >  X.  Again  using  the  fact  that  conditioning  on  less  reduces 
variance, 

Var  ( E[  I  o(wJ])  <  Var  (  e[  |  o(w0)]),  (27) 

illustrating  the  power  of  the  conditional  Monte  Carlo  approach. 

In  the  setting  of  section  2.2  specialized  to  continuous-time  Markov 
chains,  similar  arguments  give  counterparts  to  (25)  and  (27).  Since 
uniformization  increases  both  variance  and  work,  don't  use  it  there  either. 
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